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Temporal coherence is a fundamental property of macroscopic quantum systems, such as lasers 
in optics and Bose-Einstein condensates in atomic gases and it is a crucial issue for interferometry 
applications with light or matter waves. Whereas the laser is an "open" quantum system, ultracold 
atomic gases are weakly coupled to the environment and may be considered as isolated. The co- 
herence time of a condensate is then intrinsic to the system and its derivation is out of the frame 
of laser theory. Using quantum kinetic theory, we predict that the interaction with non-condensed 
modes gradually smears out the condensate phase, with a variance growing as At 2 + Bt + C at long 
times t, and we give a quantitative prediction for A, B and C. Whereas the coefficient A vanishes for 
vanishing energy fluctuations in the initial state, the coefficients B and C are remarkably insensitive 
to these fluctuations. The coefficient B describes a diffusive motion of the condensate phase that sets 
the ultimate limit to the condensate coherence time. We briefly discuss the possibility to observe 
the predicted phase spreading, also including the effect of particle losses. 

PACS numbers: 03.75.Kk, 03.75.Pp 



I. INTRODUCTION 

Bose-Einstein condensation eventually occurs in a 
bosonic system, if one reduces the temperature at a fixed 
density. It is characterized by the macroscopic occupa- 
tion of the lowest single particle energy mode and by the 
onset of long range coherence both in time and space. 
Initially predicted by Einstein for an ideal Bose gas in 
1924, it has now been observed in a wide range of phys- 
ical systems: in liquid helium [l], in ultracold atomic 
gases [1, 4] , and in a variety of condensed matter systems 
such as magnons in anti-ferromagnets Q, and exciton 
polaritons in microcavities Q . Among all these systems, 
ultracold atomic gases offer an unprecedented control on 
experimental parameters and allow very precise measure- 
ments as is custom in atomic physics. Experimental in- 
vestigation of time coherence in condensates began right 
after their achievement in the laboratory 0, [H, 0] and 
the use of condensates in atomic clocks or interferome- 
ters is currently a cutting-edge subject of investigation 
Ql ED, El, El- Therefore a crucial issue is to deter- 
mine the ultimate limits on the coherence time of these 
systems. Unlike lasers and most solid state systems in 
which condensation has been observed, ultracold atomic 
gases are weakly coupled to their environment. The in- 
trinsic coherence time of a condensate is then due to its 
interaction with the non-condensed modes in an ideally 
isolated system, which makes the problem unique and 
challenging. For the one dimensional quasi-condensate 
a theoretical treatment exists 



compared with experiment [15 



141 that was successfully 
la ]. In a true three di- 
mensional condensate, the problem was solved in [13] at 
zero temperature while until now it has been still open 
at non-zero temperature. 



As it is known since the work of Bogoliubov |l8[, the 
appropriate starting point for the description of a weakly 
interacting degenerate Bose gas is that of a weakly inter- 
acting gas of quasi-particles: the Bogoliubov excitations. 
The interactions among these quasiparticlcs shall play 
the main role in our problem. They have to be included 
in the formalism in a way that fulfills the constraint of 
energy conservation, a crucial point for an isolated sys- 
tem. A first set of works addressed the problem of phase 
coherence in condensates usin g op en-system approaches 
in analogy with the laser [l9l. |2ClL l2lj : diffusive spread- 
ing of the condensate phase was predicted. These works 
however are not to be considered as quantitative, due to 
the fact that a simplified model is used in [l9j . and due 
to an approximate expression of the condensate phase 
derivative in [2(| [2l[ ■ Moreover, lacking the constraint of 
energy conservation, these approaches neglect some long 
time correlations among Bogoliubov excitations that are 
responsible for a ballistic spreading in time of the con- 
densate phase as shown in [2 21 . |23| using many-body ap- 
proaches. Unfortunately the final prediction in [22| does 
not include the interactions among Bogoliubov modes: 
the Bogoliubov excitations then do not decorrelate in 
time, the prediction quantitatively disagrees with quan- 
tum ergodic theory [23|, and no diffusive regime for the 
condensate phase is found. Finally, the ergodic approach 
in [23[ , while giving the correct ballistic spreading of the 
phase, cannot predict a diffusive term. 

As we now explain, quantum kinetic theory allows to 
include both energy conservation and quasi-particle in- 
teractions, and gives access to both the ballistic and the 
diffusive behavior of the phase. To be as general and 
as simple as possible, we consider a homogeneous gas in 
a box of volume V with periodic boundary conditions. 
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The condensate then forms in the plane wave with wave 
vector k = 0. The total number of particles is fixed to 
N and the density is p = N/V . Let us consider the 
phase accumulated by the condensate during a time in- 
terval t: ip(t) = 6(t) — 8(0) where 8 is the condensate 
phase operator [24|. Due to the interactions with the 
Bogoliubov quasi-particles, the accumulated condensate 
phase will not be exactly the same in each realization of 
the experiment. We say that the phase fluctuates and 
spreads out in time or that the variance Vartp(t) is an 
increasing function of time. In presence of energy fluc- 
tuations in the initial state, the variance of the phase 
grows quadratically in time as already mentioned [22l.[23|. 
Quantitatively this may be seen as follows: for t — > oo, 
ip(t) ~ —p,(E)t/K where p,(E) is the chemical potential 
which depends only on the energy of the isolated system 
[23| . By linearizing p(E) around the average energy E 
for small relative energy fluctuations, one finds 
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This ballistic spreading in time of the phase is compara- 
ble to that of a group of cars traveling with different 
speeds. What happens if one reduces ideally to zero 
the energy fluctuations in the initial state ? We will 
show that the condensate phase will still spread but more 
slowly, with a diffusive motion. A precise calculation of 
the diffusion coefficient of the condensate phase in differ- 
ent experimental conditions, with or without fluctuations 
in the initial energy is the main goal of this paper. 

The paper is organized as follows. The most important 
section is the overview section |TT| there we present the 
main results of the paper that we test against classical 
field simulations, and we indicate two possible schemes to 
observe them experimentally with cold atoms. Further 
precisions and all the technical details are given in the 
subsequent sections. Starting from kinetic equations in 
section 11111 that we linearize and solve in section HV1 we 
obtain explicit results for the phase variance in section 
fVl We discuss the effect of losses in section I VII and we 
conclude in section IVlTl 



II. OVERVIEW AND MAIN RESULTS 

For a low temperature gas T <C T c the temporally 
coarse-grained derivative of the condensate phase can be 
expressed in terms of the numbers fik of quasi-particles 
of wave vector k 1231 



m being the atom mass, and Uk, Vk are the coefficients 
of the usual Bogoliubov modes: 
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As a consequence of ([2]), the variance of the condensate 
phase is determined by the correlation functions of the 
Bogoliubov quasiparticle numbers h^. Let C(t) be the 
time correlation function of the condensate phase deriva- 
tive (H: 



c(t) = (mm) - m))(0(o)) 



(4) 



By integrating formally <p(t) over time and using time 
translational invariance: 

Vaxip(t) = 2t [ C(r)dT-2 [ tC(t)(It . (5) 
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From equation ([5]) we see that two possible cases can 
occur. If C(t) is a rapidly decreasing function of r so 
that the integrals converge for t — > oo, the variance of 
the phase will grow linearly in time for long times and 
the condensate phase undergoes a diffusive motion with 
a diffusion coefficient 



D= I C(r)dT. 
Iq 



(6) 



If C(t) tends to a non zero constant value for r — > oo, 
the phase variance grows quadratically in time and the 
phase undergoes a ballistic spreading. The two differ- 
ent scenarios are illustrated in Figure Q] To describe 
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FIG. 1: Schematic view of the correlation function of the 
condensate phase derivative C(t). If C(t) tends to zero fast 
enough for t — > oo, the phase spreading is diffusive. If C(t) 
tends to a constant A, the phase spreading is ballistic. 
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where the constant term po is the ground state chem- 
ical potential of the gas and A k = j^(Uk + Vk) 2 . The 
coupling constant g for interactions between cold atoms is 
linked to the s-wave scattering length a by g = Airh 2 a/m, 



the evolution of the quasiparticles number fluctuations 
#nk(t) = ftk(t) — (nk) we write quantum kinetic equa- 
tions [26| that we linearize. Introducing the vector A of 
components A^, the vector x(t) of components Xk(t): 

afc(t) = J2 A * (*Mt)<fftk' (0)> , (7) 
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and the matrix M of linearized kinetic equations, one 
has: 

= Mx(t) . (8) 

Knowing x(t) we can calculate the phase derivative cor- 
relation function as 

C(t) = A-x(t). (9) 

On the basis of these equations we get our main result, 
that is the asymptotic expression of the variance of the 
condensate accumulated phase at long times: 

Var <p(t) ~At 2 + Bt + C for t -> oo . (10) 

In what follows we give an explicit expression for the 
coefficients A, B and C. 
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FIG. 2: Rescaled diffusion coefficient of the condensate phase 
as a function of the rescaled temperature. Full line: numer- 
ical result from the solution of pip. Dashed line: analytical 
prediction at low temperature: y = 0.3036a; 4 (see Appendix 
iDl . Dotted line: approximate prediction of linear scaling at 
high temperature: y oc x (see Appendix [Ej . 

The matrix M has a zero frequency eigenvector uq . We 
then split the correlation vector x into two components: 
x — ^UQ + X^t). The component of x along uq is constant 
in time. If it is non zero, C(t) does not decay to zero for 
t — ► oo and the phase variance will grow quadratically. 
In our general formalism we can show that 7 is linked to 
energy fluctuations in the initial state and we recover the 
result |T]) for the coefficient A. The remaining component 
X has zero mean energy. For the linear coefficient ruling 
diffusive phase spreading we find B = 2D with: 

D = -A-M~ 1 X(0), (11) 
and for the constant term 

C = -2A- M~ 2 X(0) . (12) 

Remarkably X, and thus D and C, do not depend on 
the energy fluctuations of the initial state, up to second 
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FIG. 3: (Color online) Variance of the condensate accumu- 
lated phase as a function of time for kBT/pg = 10. Black 
full lines: Var<£(t). Dashed lines: asymptotic behavior (I10|l . 
Red line (axis labels on the right): correlation function of the 
phase derivative C(t). The upper curves for Var </?(£) are ob- 
tained in presence of canonical ensemble energy fluctuations 
in the initial state. The lower curves, as well as C(t) cor- 
respond to the microcanonical ensemble where .4 = 0. In 
typical atomic condensates the healing length £ such that 
ft 2 /2m£ 2 = pg is at most in the /xm range and the unit of 
time h£ 3 1 g is at most in the ms range. 



order in the relative energy fluctuations. We find that, 
in the thermodynamic limit, the rescaled diffusion coef- 
ficient HDV/g is a universal function of ksT/pg that we 
show in Figure [2] This universal scaling was also found in 
[27I I in the frame of a classical field model. At low temper- 
ature, we have shown analytically that HDV/g scales as 
the fourth power of ksT/pg, while at high temperature 
the rescaled diffusion coefficient grows approximately lin- 
early with ksT/ pg (we expect logarithmic corrections to 
this law). As made evident by the rescaling, D is pro- 
portional to the inverse of the system volume and thus 
vanishes in the thermodynamic limit. The same property 
holds for C, and also for A in the case of canonical en- 
semble energy fluctuations. In Fig|3]for the temperature 
value ksT/pg = 10 we show the correlation function of 
the condensate phase derivative C(t), that we calculate 
integrating ijHJ) in time. On the same plot we show the 
variance of the condensate phase as a function of time 
that is obtained from ([5]). The asymptotic behavior of 
Var <p(t) from equation (|10p is reached after a transient 
time that is typically the decay time of the correlation 
function C(t). 

In the high temperature regime ksT/pg 3> 1, we were 
able to test our predictions against exact simulations 
within a classical field model. In order to perform a quan- 
titative comparison, we rephrased our kinetic theory for 
a classical field on a cubic lattice. In both the classical 
kinetic theory and the classical field simulations we intro- 
duce an energy cut-off such that the maximum energy on 
the cubic lattice is of order ksT [2a ]. We show the result 
of the comparison in Figf?] As expected the numerical 
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FIG. 4: Diffusion coefficient D c i from the classical field the- 
ory on a lattice as a function of the temperature. Crosses 
linked by a line: results from the classical version of our ki- 
netic theory. Bullets with error bars: results from the classical 
field simulations with 1000 stochastic realizations in the mi- 
crocanonical ensemble. In both curves there is a a cutoff at 
energy ksT [28j]. 

value of the diffusion coefficient D c i is different from the 
exact one given by the quantum theory and it depends 
in particular on the value of the cut-off. From the figure 
we find nevertheless a remarkable agreement between the 
classical kinetic theory and the classical field simulations. 

Our findings could have an immediate impact on 
present experiments with atomic condensates. Phase 
measurements have indeed already been successfully per- 
formed within two main schemes. 

The first scheme is out of equilibrium: starting from 
a condensate in a given internal state a, one applies two 
coherent short electromagnetic pulses separated by an 
evolution period for the condensate phase. Each pulse 
transfers a fraction of the atoms into another internal 
state b. After the second pulse one measures the num- 
ber of atoms in state b. In the original realization of this 
interferometric scheme Q, tt/2 pulses were used which 
produce a strongly out of equilibrium state of the sys- 
tem with a complex phase dynamics (29J. We propose 
to transfer only a tiny fraction of atoms in each of the 
two pulses, so that the depletion of the a condensate and 
the interactions within b atoms may be neglected. More- 
over a spatial separation of a and b 30] or a Feshbach 
resonance [HI, [3l[ may be used to suppress the a — b in- 
teractions. The ideal limiting case would be to transfer 
in 6 a single atom which could be detected in a high fi- 
nesse optical cavity [32| . Using linear response theory 
one finds that the number of atoms in b after the sec- 
ond pulse is proportional to N + $te(e lSt ° (aJ(io)ao(O))), 
where arj is the condensate operator, S is the detuning 
of the coherent pulses from the single atom a — b transi- 
tion and to is the time interval between the two pulses. 
This signal is directly dependent on Var r 3(t). Indeed 
\(al{t)a (0))\ ~ JVexp[-Var0(t)/2]. 

The second scheme uses a symmetric atomic Joseph- 



son junction [13|, |33j, in which one would cut the link 
between the two condensates by raising the potential bar- 
rier and measure the relative phase after an adjustable 
delay time. In this case an additional source of ballis- 
tic phase spreading is the partition noise proportional to 
the variance of the relative atom number. For homoge- 
neous systems with canonical ensemble energy fluctua- 
tions on both sides of the Josephson junction, the ra- 
tio between this undesired contribution and At 2 scales 
as CAr(/° a3 ) -1 ^ 2 /-^ where £jy is the number squeezing 
parameter of the Josephson junction, on the order of 
0.35 in [H, and A = A/[(pg/h) 2 (a 2 £/V)}, where £ is 
the healing length, depends only on ksT/pg 23]. For 
ksT/pg = 5 one has A ~ 150 so that, for the typical 
value (pa 3 ) 1 / 2 = 2.5 x 10~ 3 , the undesired contribution 
is smaller than At 2 [3^ ]. 



III. KINETIC EQUATIONS FOR THE 
BOGOLIUBOV EXCITATIONS 

At low temperature T -C T c we assume that the state 
of the gas can be approximated by a statistical mixture 
of eigenstates of the Bogoliubov Hamiltonian -ffeog 

H Bog = E + ^2 e k h^ , (13) 

where Eq is the energy of the ground state. The eigen- 
states of i?Bo g are Fock states |{^k}) with well defined 
numbers of Bogoliubov quasiparticles. Whereas expec- 
tation values of stationary quantities are expected to 
be well approximated by Bogoliubov theory, this is no 
longer the case for two-time correlation functions. This 
is physically quite clear for the correlation function of the 
Bogoliubov mode occupation numbers h-^: whereas they 
never decorrelate at the Bogoliubov level of the theory 
(they are conserved quantities of ^Bog), they will experi- 
ence some decorrelation for the full Hamiltonian dynam- 
ics because of the interactions among Bogoliubov quasi- 
particles, that are at the origin of the Beliaev-Landau 
processes. 

For a given initial state of the system characterized by 
the occupation numbers {n^} = {fik(0)} the time evo- 
lution, beyond Bogoliubov approximation, of the mean 
mode occupation numbers 

«q(t) = ({nk(0)}|n q (t)|{n k (0)}) (14) 

can be described in terms of quantum kinetic equations 
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of the form 2C 



9 P 
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In ()15|) we have introduced in (|15p the coupling ampli- 
tudes among the Bogoliubov modes: 



Al, k , = U q U k U k , +V q V k V k , 

+ {U q + V q )(V k U k ,+U k V k ,). 



(16) 



Kinetic equations (fT5|) describe Landau and Beliaev pro- 
cesses in which the mode of wave vector q scatters an 
excitation of wave vector k giving rise to an excitation 
of wave vector k' (Landau damping) , the mode of wave 
vector q decays into an excitation of wave vector k and 
an excitation of wave vector k' (Beliaev damping) , and 
inverse processes. In each process the final modes have 
to satisfy energy and momentum conservation. Energy 
conservation is ensured by the delta distributions in (|15|) 
where e k is the Bogoliubov energy of the quasiparticle of 
wave vector k, 
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To calculate the correlation function C(t), equations 
(|15|) can be linearized for small deviations [35| . and lin- 
ear equations for the correlation functions Xq(t) can be 
obtained: 



x = M x . 



(18) 



To obtain i q from n q , we connect expectation values in 
an initially considered Fock state to expectation values in 
the system state by an additional average. More details 
on the derivation of (|18[) . as well as the explicit form 
of the equations, which arc in fact integral equations, 
are given in appendix |XJ In particular, the matrix M 
depends on the Bose occupation numbers 
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The set of h q constitutes a stationary solution of (fT"5|) . 
with a temperature T such that the mean energy of this 
solution is equal to the mean energy of the system. The 
classical version of kinetic equations that we used to test 
our results against classical field simulations (that are 
exact within the classical field model) are reported in 
appendix [B] 



IV. SOLUTION OF THE LINEARIZED 
EQUATIONS 

The matrix M is real and not symmetric. It has right 
and left eigenvectors Mu\ = m\u\, t v\M = ni\ t v\ sat- 
isfying v\ ■ uy = 5\\: . Due to the fact that the system 
is isolated during its evolution, M has a pair of adjoint 
left and right eigenvectors with zero eigenvalue [36j. In- 
deed for any fluctuation 5n, introducing the vector e of 
components e k , one has 



y^e fc n k 

k 



constant — > } ekSn-^ = — > *e MSn = . 

(20) 

Let us denote uo the right eigenvector of M with eigen- 
value and vq the corresponding left eigenvector. One 
has from (f2"Dl) v = e. On the other hand one can show 
that uo = a [37j with : 



a k 



£kn k {n k + 1) 



(21) 



It is useful to split the correlation vector x into a com- 
ponent parallel to a and a zero-energy component, that 
is a component orthogonal to the vector e: 



"fa - 



X. 



(22) 



For our normalization of a one simply has 7 = e-x. From 
equations (fT5f and (|22|) we then obtain 



7 
X 





MX. 



(23) 
(24) 



Under the assumption that A- X{t) = 0(t~( 2+ ^) with 
v > for t — ► 00, we obtain from ([5]) the asymptotic 
expression for the condensate phase variance: 



Yai(p{t) 
with 



A 
B 



■Bt + C + o(l) for i->oo (25) 



a ■ jd 



drA ■ X{t) 



poo 

C = -2 dr tA ■ X(t) , 
Jo 



(26) 
(27) 

(28) 



As explained below, in the paragraph "The correlation 
function C(t)" of section [Vj we have some reason to 
believe that A ■ X(t) scales as t~( 2+ ") for large r with 
v=\. 



V. RESULTS FOR THE PHASE VARIANCE 

State of the system and quantum averages: In the 

general case, we assume that the state of the system is 
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a statistical mixture of microcanonical states. For any 
operator O one then has 



(O) = / dEP(E)(6) mc (E), 



(29) 



where (. . .} mc (£) is the microcanonical expectation value 
for a system energy E. Furthermore we make the hy- 
pothesis that the relative width of the energy distribution 
P{E) is small. Formally, in the thermodynamic limit we 



E 



for N — > oo . 



(30) 



Besides microcanonical averages (0) mc (E), we introduce 
canonical averages (0) cain (T) where the temperature T 
is chosen such that (H-Q 0g ) can (T) = (#Bo g ) = E. Useful 
relations among the quantum averages in the different 
ensembles are derived in Appendix [Cl 
Quadratic term: First we calculate the quadratic term 
A of the condensate phase variance given in (|26|) . We 
introduce the "chemical potential" operator 



fi = + H ^ ^k/lk 
k#0 



(31) 



so that —fi/h = (p according to equation ([2|). The con- 
stant 7 appearing in (|26[) can then be expressed as 



7 = €■ x (0) = < (tf B og - E) fi)/h (32) 
so that, using (|29| . 



Finally 



A = 



V&rE. 



(38) 



We then recover, by a different method and in a more 
general case, the main result of [23[ for super diffusive 
phase spreading when energy fluctuations are present in 
the initial state of the gas. 

Linear term: The linear term B in (j27|) represents a 
diffusion of the condensate phase with a diffusion coeffi- 
cient D = B/2. Integrating equation from zero to 
infinity and assuming X(oo) = 0, we obtain 



D 



-1-M- 1 X(0) 



(39) 



where the inverse of the matrix M has to be understood 
in a complementary subspace to the kernel of matrix M, 
that is in the subspace of vectors x satisfying e- x = 0. 
We can then write 



D = -{PA) ■ M _1 X(0) 



(40) 



where the matrix P' projects onto this subspace in a 
parallel direction to a. This corresponds to a matrix P 
given by 

-Pk.k' = <Sk,k' - ekak' ■ (41) 
As a consequence, one simply has 

X(0) = ptf (0) = f(0) - a (e ■ f(0)) . (42) 



7- dEP(E)(E-E)(fi) mc (E)/h. 



(33) 



We now expand the function (fl} mc (E) around its value 
for the average energy: 



(A)mc(£) = (A)mc(£) + {E ~ E) 



dE 



-(E) 



(34) 



Inserting the expansion (|34j) in (|33|) one gets to leading 
order in the energy fluctuations: 



d(fi) mc - V&tE 
^^E~ iE) — 



(35) 



Using equation (|C4|) of Appendix [Cl for O — fi, we finally 
obtain 



7 ~ i^blVari?. 



(36) 



According to (|2"6")l we also need the value of A ■ a that we 
can rewrite using (|C6|) . (|C7[) as 



d 77- 



h 4r E 



(37) 
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FIG. 5: Constant C as a function of the rescaled temperature. 
Full line: numerical result from the solution of (|52p . Dashed 
line: analytical prediction at low temperature: y — 0.2033a; -1 
(see Appendix |Dj. Dotted line naive prediction for the high 
temperature scaling: y oc x 1 ^ 2 (see Appendix |E)| , 

We show here that D does not depend on the width of 
the energy distribution P(E) of the initial state. To this 
end it is sufficient to show that the same property holds 



7 



for X(0). We apply the relation (|C10[) to and fi\ji^ 
to obtain after some calculations 

(<5n k <5n k ') = 8^n k {l + n k ) 

+ (r, - l)fc B T 2 ^^"^ + ...(43) 

dT E 

where the dots indicate terms giving higher order con- 
tributions in the thermodynamic limit that will be ne- 
glected. Here r] is the ratio of the variance of the system 
energy to the energy variance in the canonical ensemble, 
i] = Var E /\ai ca _ n E. Eq. ([43|) shows that x(0) and hence 
X(0) are affine functions of r\. X(0) can then be deter- 
mined from its values in r\ = (microcanonical ensemble) 
and 77 = 1 (canonical ensemble): 

X(0) = ryX can (0) + (1 - T])X mc {0) . (44) 

On the other hand one can show explicitly for a large 
system that X can (0) = X mc (0) [38| . As a consequence 

X(0)=X mc (0) (45) 

does not depend on rj. Note that this relation extends 
to all positive times, X(t) = X lnc (t), since the matrix M 
does not depend on the energy fluctuations. 

The expression of X mc (0) has been derived in [27| . In- 
troducing the covariance matrix of Bogoliubov occupa- 
tion numbers 

QS,(f)Hfc k (t)%(0)), (46) 

one has in the microcanonical ensemble 

^mc(O) =Q mc (i = 0)l. (47) 

As we showed in 27], for a large system, the t = 
covariance matrix in the microcanonical ensemble can be 
obtained by the one in the canonical ensemble by projec- 
tion: 

Q mc (t = 0) ~ ptQ can (yj = 0)P , (48) 

where Q can is the covariance matrix in the canonical en- 
semble, that can be calculated using Wick's theorem 

Qk£(* = 0) = n fc (n fc + l)$i t>k ,. (49) 

Using (|40|) we can then calculate the diffusion coefficient 
D already discussed in the paper and shown in FigfS] 
Some details about the low temperature and high tem- 
perature limits of D are given in appendix [D] and in ap- 
pendix |E] respectively. In particular we find at low tem- 
perature 

^~ C1 (MY for M^o (50) 
a \ pa ) pa 

The constant ci = 0.3036 is calculated numerically. 




» L , 1 , 1 , 1 , 1 , 1-4x10' 

'0 200 400 600 800 1000 

t/(nlf/g) 




t/(fl^/g) 

FIG. 6: (Color online) Top: For a system prepared in the 
microcanonical ensemble, variance of the condensate accu- 
mulated phase as a function of time for ksT/pg = 0.2. Black 
full line: Var<£ obtained from (0). Dashed line: asymptotic 
behavior (|25[1 . Red line: correlation function C(t) defined 
in Q. Bottom: Full red line: The same correlation func- 
tion C(t) in log- log scale. Dotted line: Exponential function 
f(t) — C(0) exp(— t/r c ), where r c is defined in (|54[1 . Dashed 
dotted line: law y oc x~ 3 predicted by the Gaussian model of 
[13 • £ is the healing length: h 2 /2m£, 2 = pg. 

The constant term: We now come to the constant term 
C defined in (28]), By integrating formally (d/dt)(tX) 
between zero and infinity and by using (|24p. we obtain 

/>oo />oo 

0= / dtX(t)+M dttX(t) (51) 
Jq Jo 

and finally 

C = -2{PA) ■ M- 2 X{0) . (52) 

We show in Fig[5]the constant C obtained from (|52|) as a 
function of temperature. At low temperature we get 

CV fkeTy 1 k B T 

-pr ~ c 2 for > (53) 

£ d V pg J pg 
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The constant C2 = —0.2033 is calculated numerically. 
Note that, contrarily to the coefficients A and B, the 
coefficient C does not tend to zero for T — > 0, on the con- 
trary it diverges. However, the typical decay time r c of 
the correlation function C(t) also diverges in this limit, 
as we shall see in what follows. 

The correlation function C(t): The phase derivative 
correlation function C(t) was defined in (TJ|. Restricting 
for simplicity to the system being prepared in the mi- 
crocanonical ensemble (as we have seen, in the general 
case, C(t) deviates from the microcanonical value C mc (i) 
by an additive constant), we show in FigJBJ-Top the func- 
tion C mc (t) in the low temperature case ksT/pg = 0.2. 
C mc (t) is obtained by ©, integrating equation (fTgj) in 
time by Euler's method. Correspondingly, we calcu- 
late the variance of the condensate accumulated phase 
as a function of time from ([5]) and we compare it to its 
asymptotic behavior (|2"5|) . On the same figure, see FigO 
Bottom, we show C mc (t) in log-log scale to point out sig- 
nificant deviations from the exponential behavior: C(t) 
rather decays as a power law; the Gaussian model of 27 1 
at large times gives C mc (t) oc 1/t 3 which we also plot in 
the figure for comparison. 

Characteristic time to reach the asymptotic 
regime: The asymptotic regime for the phase variance 
is reached after a transient that is the typical decay time 
of the correlation function C(t). An estimation of this 
time is 



D 



Cmc(0) 



(54) 



This is only an estimation since, as we have seen, C mc (i) 
is not an exponential function oc exp(— t/r c ). A plot of 
r c as a function of temperature is shown in FigJT] At low 
temperature 



9Tc 



C3 



k B T 

pg 



for 



k B T 

pg 



(55) 



The constant C3 = 0.05472 is calculated numerically. 

In table Q] we give the numerical values of the rele- 
vant parameters for 10 reduced temperatures in the range 
0.1 - 100. 




k B T/pg 

FIG. 7: Typical decay time of the correlation function C(t), 
after which the asymptotic behavior of the phase is ob- 
served. Full line: numerical result from the solution of 
(|54[) . Dashed line: analytical prediction at low temperature: 
y — 0.05472a; -5 (see Appendix IT)]) . Dotted line: naive pre- 
diction for the high temperature scaling: y oc a; -1 / 2 (see Ap- 
pendix [E}. The healing length £ is such that ft 2 /2m£ 2 = pg. 



TABLE I: Numerical values of the relevant quantities. £ is 
the healing length: h 2 /2m£ 2 = pg. A is given for energy 
fluctuations of the canonical ensemble, and (7(0) for the mi- 
crocanonical ensemble. 



k B T 


hDV 

a 


cv 
W 


C mc (0)Vft 2 £ 3 


9T C 

hF 




pg 




(pg)a^/ 2 


0.1 


2.130 x 10~ 5 


-2.227 


1.784 x 10~ 9 


11940 


0.02397 


0.2 


2.142 x 10 -4 


-1.426 


2.046 x 10~ 7 


1046 


0.1092 


0.5 


3.163 x 10~ 3 


-1.286 


3.105 x 10~ 5 


101.9 


0.6037 


1 


1.911 x 10 -2 


-1.726 


6.337 x 10" 4 


30.16 


1.7557 


2 


9.626 x 10~ 2 


-2.886 


7.939 x 10~ 5 


12.12 


4.3682 


5 


0.638 


-6.691 


0.134 


4.746 


12.276 


10 


2.280 


-12.95 


0.880 


2.590 


24.542 


20 


7.323 


-24.48 


4.971 


1.473 


46.598 


50 


30.14 


-53.35 


42.06 


0.716 


103.10 


100 


81.60 


-91.59 


195.8 


0.417 


182.94 



VI. INFLUENCE OF PARTICLE LOSSES ON 
THE SUPER-DIFFUSIVE PHASE SPREADING 

For an isolated system with energy fluctuations in the 
initial state, we have seen that the correlation function 
C (t) of the condensate phase derivative does not vanish 
at long times and the condensate phase spreading is super 
diffusive. In presence of particle losses, unavoidable in 
real experiments, the system is not isolated and the total 
energy is not conserved so that one may wonder whether 
the super diffusive term is still present. We show in this 
section that this is indeed the case, in a regime where the 
fraction of particles lost during the decay time r c of the 



correlation function C(t) is small, a condition satisfied in 
typical experimental conditions. 

We first perform a classical field simulation with one 
body losses of rate constant T: during the infinitesimal 
time interval dt, a quantum jump may occur with a prob- 
ability TNdt where N is the number of particles just 
before the jump. If the jump occurs, a particle is lost 
which corresponds in the classical field model to a renor- 
malization of the field ip(r) -► [{N - l)/N] x ' 2 ii{r). In 
between jumps the field evolves with the usual non-linear 
Schrodinger equation: 

h 2 

ihd t ^ = -—A^ + g\i;\ 2 il>. (56) 
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This results from the interpretation of the classical field 
in terms of an Hartree-Fock Ansatz for the quantum sys- 
tem state as detailed in Appendices iFl and [Gl The result 
for the condensate accumulated phase standard deviation 
as a function of time is shown in Fig. [8j in the absence 
(dashed line) and in presence (solid line) of losses. It is 
apparent that, for the parameters taken in this figure, 
the spreading of the phase up to a standard deviation of 
order unity is only weakly affected by the particle losses. 
We also find that the phase spreading is in fact acceler- 
ated by the losses and becomes effectively super-ballistic. 
As we now show, this is due to the fact that the losses 
introduce particle number fluctuations that grow in time. 




5000 



t/(*£7g) 



FIG. 8: Condensate accumulated phase standard devia- 
tion as a function of time with and without one body 
losses in a classical field model. Solid line: simulation for 
r = 1.555 x 10~ 5 g/(h(, 3 ). Dashed line: simulation with- 
out losses. Black discs: lossy ergodic model (see text) for 
r = 1.555 x 10~ 5 g/(h£ 3 ). Circles: prediction of the ergodic 
theory (no losses). The initial atom number is iV(0) = 4x 10 , 
p(0)g = 1798.47h 2 /(mV 2/3 ), k B T/p(0)g = 2.95. A spatial 
box of sizes Li, L2, L3 and volume V = L1L2L3 is used with 
periodic boundary conditions. The squared box sizes are in 
the ratio y/2 : (1 + VE)/2 : \/3- Note that here, contrarily to 
previous figures, the variance is directly given and was not di- 
vided by the factor ( 3 /V (here £ 3 /V ~ 4.64x 1CT 6 ). For a typ- 
ical atomic density of p(0) = 1.2 x 10 20 atoms/m 3 , taking the 
87 Rb mass and scattering length a = 5.3 nm, our parameters 
correspond to 1/T = 20s, p(0)g/(2-irh) ~ 950 Hz, T ~ 0.14juK 
or T ~ 0.3T C , and the temporal unit H£ 3 /g — 0.31ms. A num- 
ber of 1200 realizations is used in each simulation, and the en- 
ergy cut-off corresponds to a maximal Bogoliubov eigenenergy 
equal to ksT. The variance of the total energy in the initial 
state is 1.5 x 10 11 h 4 / (mF 2 ' 3 ) 2 , resulting from sampling the 
canonical ensemble in the Bogoliubov approximation. This 
value is larger than the one predicted by the Bogoliubov the- 
ory by a factor 1.3 due to non negligible interactions among 
the Bogoliubov modes. A lossless relaxation phase of a dura- 
tion 500ft£ 3 /<? is performed after the sampling. 

In order to understand the numerical results we use a 
heuristic extension of the ergodic model in presence of 
losses. In the model there are two dynamical variables: 
the total energy and the total number of particles. We 



assume that in between two loss events the condensate 
phase evolves according to 



9(t) = 



(57) 



where /i mc is the chemical potential in the microcanon- 
ical ensemble of energy E for a system with N parti- 
cles. When a loss event occurs, N is obviously changed 
into N — 1. For the energy change one has to con- 
sider separately the kinetic and the interaction energies: 
i?km is a quadratic function of ip and is changed into 
[(N— 1) /N] i?kin • The interaction energy is a quartic func- 
tion of tp and is changed into [(N — 1)/A] 2 i?i n t. When a 
jump occurs we then take 



E' 



N' = 



N - 1 

N 

N - 



N 
N -1 



(E kin ) mc (E,N) 



(E- mt ) mc (E, N) 



(58) 
(59) 



where the prime indicates the quantities after the jump 
and where (-Ekin)mc and (£int)mc are the mean kinetic 
and interaction energies in the microcanonical ensemble 
[3^ | . To calculate the microcanonical averages and the 
chemical potential, we rely on Bogoliubov theory. In the 
classical field model: 



^nc(E,N) 

(E int ) mc (E,N) 
(E kin ) mc (E, N) 



gN g E-E 
V V M 



E 



Eo 



gN E-E a 



V M 

E — (£'int)mc 



E 

k#0 



h 2 k 2 

2m 

h 2 k 2 

2m 



2gN 
V 

2gN 



60) 

)ei) 



(62) 



where M. is the number of Bogoliubov modes and Eq = 
gN 2 /2V is the ground state energy. 

We have performed a Monte Carlo simulation of this 
model. The initial energy is obtained sampling a Gaus- 
sian distribution with a mean energy given by Bogoliubov 
theory and with the same variance as in the classical field 
simulations (see caption of Fig[8]). The results for the 
condensate phase variance (symbols) are compared with 
the classical field simulation with and without losses in 
FiglU A good agreement is found. 

To go further, we analytically solve this model to first 
order in the loss rate constant T. As detailed in appendix 
IFl we obtain the simple result: 



Vaxtup(t) ~ (Var^i)t 2 



(( M <5 M )-( M )(^))r7Vt 3 
1 



2 TNt 3 . 



(63) 



Here N is the initial atom number, /z is the initial chem- 
ical potential, and Sfj, = fi! — /i is its change after the 
first loss event. After an explicit calculation, in the limit 
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ksT pg, this reduces to 



Var <p(t) 



dp r _ 



dE 



(E,N)VaxE 



h 2 



-FNt 3 (—Y 
3 \hv) 

\hv) \ N 



(64) 



where (SN) /N -C 1 is the fraction of non condensed par- 
ticles. The first term in the right hand side of ([64]) is 
the classical field version of the result ((T|) without losses. 
The third term is negligible as compared to the sec- 
ond one in the present regime of a small non condensed 
fraction. The second term, independent of the temper- 
ature, is the result that one would have at zero tem- 
perature in presence of losses at short times (Tt -C 1). 
This term has a simple physical interpretation: in pres- 
ence of fluctuations in the initial number of particles 
for a lossless pure condensate, the condensate accumu- 
lated phase grows quadratically in time with a variance 
(dp,/dN) 2 t 2 Va,r N/H 2 , where p = gN/V. For a lossy pure 
condensate with initially exactly N particles, Var TV ~ 
FNt so that one indeed expects Var tp(t) cx (g/HV) 2 TNt 3 . 
Actually, at T = it is possible to calculate exactly 
Var tp(t) in presence of losses (see appendix UJ: 



[Vax^)] T=0 = y^) 2 v[l 



2Tte 



-rt 



-2rti 



(65) 

This zero temperature result even extends to the quan- 
tum case for a pure condensate, see Appendix [Ul so that 
one may hope that the form of the classical field result 
extends to the quantum reality. To be complete 
we also give the exact value of the correlation function 
(cio(f)ao(O)) in the quantum case for a pure condensate 
with initially TV particles and subject to one-body losses: 



(oJ(t)oo(O)) 



T=0 



-rt/2 N 



r 

A 



n N-l 



-At\ 



(66) 

with A = r — ^7. This can be obtained by applying 
the quantum regression theorem using (jGlj) expressed in 
the Fock basis. The same result (|66|) can be obtained 
using the exact result for a two mode model Eq.(125) of 
[4~fl | , and assuming that the second mode of infinitesimal 
population experiences no interactions and no particle 
losses. 



VII. CONCLUSIONS 

In conclusion we have presented a full quantitative 
quantum solution to the long standing problem of the 
decoherence of a condensate due to its interactions with 
quasi-particles in the non-condensed modes: the growth 
of the variance of the condensate accumulated phase in- 
volves in general both a quadratic term and a linear term 



in time, with coefficients that we have determined within 
a single theoretical frame, quantum kinetics. As we have 
discussed, our findings may be directly tested with state- 
of-the-art technology, and they may stimulate systematic 
experimental investigation of this problem, both funda- 
mental and crucial for future applications of condensates 
in matter wave interferometry. 
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APPENDIX A: EQUATIONS FOR x<,(t) 



We detail here the derivation of equation ([18]) for the 
correlation functions Xq(t). We assume that (i) the den- 
sity matrix p of the gas is a statistical mixture of eigen- 
states of the Bogoliubov Hamiltonian H-Q 0g given by dE 



£7>({n k })|{ra k }><{n k }|, 



(Al) 



and 
the 



(ii) for 
evolution 



a given initial Fock state |{^k})j 
of the expectation values n q (t) = 
({n-k(0)}|n q (t)|{nk(0)}) are given by the kinetic equa- 
tions JT5]). We then have 

«(t)n q2 (0)) = ^(W(0)}K 2 (0) 

{«k(0)} 

x ({n k (0)}|n qi (t)|{n k (0)}>, (A2) 



and 



dt 



(£n qi (i)5n q2 (0)) = ]T 7>({n k (0)})fc q2 (0) 



{"k(0)> 

x X^ qi ,q'tov(*). (A3) 



where the matrix M is obtained by linearization of equa- 
tions (p~5|) . We have introduced 



Sn^t) = ({n u (0)}|n q (t)|{n k (0)}) - {h q 



(A4) 



where we recall that (. . .) is the expectation value in the 
state of the system. By multiplying (|A3|) by A q2 , sum- 
ming over q2, and approximating (h q ) with n q of (|19p . 
which is justified in the present regime of large system 
size and weak relative energy fluctuations, we obtain (|18[) . 

Using the rotational invariance of function of k 

and the delta of conservation of energy we can explicitly 
integrate over the angular variables and we obtain the 
simple integral equations that we now detail. We intro- 
duce dimensionless quantities Q. Momenta are rescaled 
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by the inverse of the healing length £ = (ti 2 /2mpg) 1 / 2 , 
energies are rescaled by the Gross-Pitaevskii chemical po- 
tential pg, and rates are expressed in units of gj (27r 2 £ 3 /i): 



( h2 V /2 p 



\ 2mpg t 

h 

99 
2tt 2 £ 3 ^ i 

5 



e 9 - ^ = [<f(<f + 2)] 1 / 2 



r, = 



(A5) 
(A6) 

(A7) 



one has 



HDV 

9 
CV 

Cmc(^) 



Pdk(PA) k (ti- 1 X(0)) k (A18) 



-2 / k 2 dk(PA) k (M-'X(0)) k (A19) 
Jo 



APPENDIX B: CASE OF A CLASSICAL FIELD 



As a consequence, the mean occupation number n q is a 
function of q and of the ratio ksT / pg only, and the mode 
amplitudes U q , V q are functions of q only. Expressing the 
time in reduced units, we then have 

x q (t)=-t q x q (t)+I. (A8) 

The integral I is: 

I f°° „- / ik < \ 2 fc(e fc + t q ){n k > - n q ) 



2- 



dk (A k k \ q )' 



q(k' 2 + 1) 



-x k (t) 



q(k" 2 + 1) 
/o V M i q{k" 2 + l) 



with 



a? fc (t)j(A9) 

(A10) 
(All) 



A: 



//2 



l + (e fc + e ? ) 2 -l 
1 + (4 - e q ) 2 - 1 . 



The damping rate T q is the sum of the Beliaev and Lan- 
dau damping rates already given in (27j : 



f = f L 4- f B 

1 H Q 



with 



l' L r +oc - r k , ^ k{l k + l q )(n k -n k >) 



g 
2tt 



dk 



{<,)' 



q(k' 2 + 1) 



(A12) 



(A13) 



and 



q - 1 dk 



( A k,k") 



2 k{l q - e k )(l + n k + n k ») 



q(k" 2 + 1) 



Introducing 



2tt 2 £ 3 /j 

M = — M 

9 

A = — A 
9 

X(t) = —X(t) 
9 



(A14) 

(A15) 
(A16) 
(A17) 



We consider a discrete model for a classical field tp( r ) m 
three dimensions. The lattice spacing is I along the three 
directions of space. We enclose the field in a spatial box 
of volume V with periodic boundary conditions. Then 
the field can be expanded over the plane waves 



(Bl) 



where k is restricted to the first Brillouin zone, k £ T> = 
[—ir/l,ir/l[ 3 . The lattice spacing corresponds to an en- 
ergy cut-off such that the highest Bogoliubov energy on 
the lattice is efe max — ksT. 

The classical limit in the kinetic equations is obtained 
by taking: n k + 1 ~ h k — > n c k l — ksT/e k in the equation 
(|A8p for Xq. In the units already introduced in Appendix 
[A] one then has: 



± q (t) 

We have introduced 

Id. = 
2 



-f q 1 x q (f)+/ cl . 



(B2) 



,3 k (Al k )\nt - nf)S(e k + e k , - e q )x k (t) 
+ J d 3 k (A k Q\ni, - nf)S(e k + e q - e k „)x k (t) 
+ I d 3 k (A h k ,, q ) 2 (nf, + n c q x )8{e k , + l q - e k )x k {t) . 



(B3) 

The integrals are restricted to the domain T> = 
{-ttP/1,t:P/1[ 3 and 



k' = q-k+yii, n e Z 3 
2tt 



(B4) 



k" = qlkljm, meZ 3 , (B5) 

where m and n are such that k',k" e V. Indeed the 
presence of the lattice implies the existence of unphysical 
Umklapp processes, such that n / or m / (see [27|). 
that we include in the classical kinetic theory. 

The damping rate in the classical field model r q is the 
sum of Beliaev and Landau damping rates r q = r q ' B + 
f q '' L with: 



= / d A k 



) 2 ("* + n «) S & + ^ ( B6 ) 
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f = 2 / d 3 k 
Jv 



(■AC) ~ nt„)5{e k + e q - l k n 



From the kinetic equations in the classical model, 
x(t) = M c \x(t), one has the classical diffusion coefficient 
in the form: 



hD cl v 

9 



(B8) 



Paradoxically the lattice with the relatively low energy 
cut-off breaks the spherical symmetry of the problem 
making the numerical solution heavier than in the quan- 
tum case. 

The classical field simulations were performed as in 
[2?J on a lattice with a few percent anisotropy, except 
for the free dispersion relation of the matter wave on the 
grid: here the usual parabolic dispersion relation Ek = 
h 2 k 2 /2m was used. 



APPENDIX C: STATE OF THE SYSTEM AND 
QUANTUM AVERAGES 

In this appendix we establish some useful relations 
among different averages. In particular we wish to ex- 
press the expectation value of O defined in equation (j2"!)|) 
in terms of canonical averages where the temperature T 
is chosen such that (-HBog)can(T) = (#Bog) = E. First 
of all we expand the function (0) mc (E) around its value 
for the average energy: 



(0) mc (E) = (0) m c(E) + (E — E) 



d{d)„ 
dE 



-(E) 



(E) 



(CI) 



We then take the average of (|C1[) over the energy distri- 
bution P(E) and obtain 



(O) = (0} mc (E) + - 



Id 2 (6), 



2 dE 2 



: Var(ff B og) 



(C2) 



The coefficient in front of Var(#Bog) in the second term 
in (|C2[) appears in a first order correction, it can thus be 
calculated to lowest order in the inverse system size. By 
writing explicitly 



(0) mc (E(T)) ~ (0) can (T) 



(C3) 



and taking the derivative of this relation with respect to 
the temperature T, we obtain 



dE 



-(E(T)) 



dE(T) 
dT 



(C4) 



and 
d 2 (0>, 



dE 2 



mn^'^t ^tfA . (C5) 



dE(T) dT \ dE(T) 
dT \ dT 



On the other hand we know that 
dnk _ 1 

It ~ k B r 2 



eknk(l + n k ) (C6) 



Var can (-ffBo g ) = ^ e 2 k n k (l + n k ) 
,dE 



knT< 



dT 



We then obtain the equation 



k B T 2 d 



dE(T) 
dT 



with 



V 



Var(g B og) 
Var can (J?Bo g ) 



(C7) 



, (C8) 



(C9) 



In the particular case in which the average (O) is taken 
in the canonical ensemble, rj = 1 and we recover equa- 
tion (B7) of [23| . If we now eliminate the microcanonical 
average in (|C8|) in favor of the canonical one, we obtain 
the final formula 



(6) ~ (0) can (T) + 



k B T 2 d 
2 dT 



dT 



(0) can (T) 



dE(T) 
dT 




APPENDIX D: LOW TEMPERATURE 
EXPANSION 



Let us consider the limit 
k B T = 
99 



e < 1. 



(Dl) 



In this case the occupation numbers n„ are exponentially 



small unless e q < e: indeed 



e /3e, _ 1 e e q /e 



1 



(D2) 



We can then restrict to low energies and low momenta 
where the spectrum is linear 



e ~ V2q for 



0. 



We thus introduce 



- q 
q = - 



e y/2k B T 



(D3) 



(D4) 



that is a dimensionless momentum of order unity for typ- 
ical Bogoliubov mode energies of order k B T. To obtain 
an expansion for £ <C 1, we then expand the relevant di- 
mensionless quantities in powers of q which is of order 
e: 



e q = V2qe + ^-q 3 e 3 + 0(e 5 ) (D5) 



(U g + V q ) 2 



v2 - ^.fe 3 + 0(e 5 ). (D6) 



-qe 
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For a general function F(f3e q ), as for example a function 



APPENDIX E: HIGH TEMPERATURE 



F(e q /e) = F(V2q) + ^f F '(V2q)s 2 + 0( £ 4 ) , (D7) 



and for the coefficients A in equation (|A9[) 

At = ^qk{ q+ ~k)E^ + 0{e^) (D8) 
' qk(q - k) e 3/2 + 0(£ 5/2 ) . (D9) 



A q 



3 
3 

2^4 



On can then write the low temperature version of equa- 
tions fA9]) . (TA1~3|) and (|A14l) : 

Q /*00 

J — e 5 — / dkk 2 (k + q) 2 (n k+q - n q )x k 
4 Jo 



+ e° 



;9tT 



(ifcfc (fc - q) {n k _ q +n q + l)x k 



+ e 5 ^ P dkk 2 {q-k) 2 {n q ^ k -n q )x k (D10) 
4 Jo 



Tl 



:9tT 



dkk 2 (k + q) 2 (n k -n k+q ) (Dll) 



o 



f B ~ e 5 ^ rdfcfc 2 (g-fc) 2 (n fe +n g _ fe + l),(D12) 

° JO 

where ~ stands for mathematical equivalence in the limit 
e — ► (/ ~ g if //<7 — * 1). In order to obtain the scal- 
ing with e of the diffusion coefficient D and of the other 
quantities, we expand PA and A(0): 



(D13) 



(Pl) q = ^ [^-g 3 ]e 3 + 0(e 5 ) 
X(0)„ = ^j-F(V2q) [qK - f] e 3 + 0(e 5 ) (D14) 



with 



F(j3e q ) = n q (n q + l) 

= J^dkk s F(y/2k) 
J °° d~k~k 4 F(V2~k) 

We then conclude that for e — > 

' Ci £ 



5 
£ 3 

9 2 



C2 S 

c 3 e" 

Q 

C4 e 



(D15) 
(D16) 

(D17) 

(D18) 
(D19) 

(D20) 



In (|D17|) - (|D20[) a factor e 3 comes from dkk 2 in the Ja- 
cobian. The numerical coefficients c\ to C4 can be calcu- 
lated numerically using the expanded expressions (jDlOP - 
(|D14p , or using the original expressions and extrapolating 
the result for k B T/pg — > 0. 



Let us now consider the high temperature limit 
k B T 



pg 



> 1. 



(El) 



A naive approach then consists in replacing the disper- 
sion relation of quasiparticles by the free particle one 



h = Vf(Q 2 + 2) - f for g^co, (E2) 



and introduce the rescaled dimensionless momentum 

q 



<1 : 



k B T 

pg 



so that q ~ e q /k B T. In this limit U q 



A lk' ^ 1, and 



1, v q 



e 0e q _ I 



1 



(E3) 

0, 

(E4) 



To lowest order, the integral I and the rate f are 

oc (k B T/ pg) 1 / 2 while X(0) and PA do not depend on 
k B T/pg. Similarly to the low temperature limit one 
could then deduce the high temperature scaling of the 
relevant quantities. However, in this naive approach in- 
frared logarithmic divergences appear in the integrals. 
By general arguments we then expect logarithmic cor- 
rections to the deduced scaling for k B T/pg — > 00. We 
can then only say that roughly 



HDV 
9 

CV 
~W 

gr c 

g 2 



k B T 

pg 

k B T 

pg 

k B T 

pg 

k B T 

pg 



1/2 

-1/2 

3/2 



(E5) 
(E6) 
(E7) 
(E8) 



where in (jE5] ) -([E8 l) a factor (k B T/pg) 3 / 2 comes from the 
Jacobian. A consequence of (|E5|) would be that at high 
temperature the diffusion coefficient is independent on g. 

This is compatible with the naive expectation that at 
high temperature the damping rate is proportional to 
the scattering cross section a — 87ra 2 oc g 2 , with a pro- 
portionality factor independent of a (as is the case for 
a classical gas where Y ~ pav where p is the density 
and v the average velocity). In this naive expectation, 
r _1 oc g~ 2 compensates the contribution of A 2 cx g 2 , 
and D is independent of g. This is actually too naive 
and neglects logarithmic corrections. For example, for 
the Landau damping rate, we were able to show that, in 
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the limit of a vanishing pg for fixed temperature T and 
momentum q: 



q P9 



-o 



(E9) 



APPENDIX F: SOLUTION OF THE LOSSY 
CLASSICAL FIELD ERGODIC MODEL 



We first derive (|65|) , by solving the lossy ergodic model 
exactly at zero temperature, and then derive (|63|) for 
T 7^ 0, by solving the lossy ergodic model to first order 
in the loss rate constant T. After an explicit calculation 
we then obtain (|64| . 

To this end, we start with the fully quantum model, 
defined by a Lindblad form master equation including 
one body losses, and we use the formulation given in 
[40| of the Monte Carlo Wavefunction method for the 
expectation value of an observable O: 



(6){t) = J2 I dh...dt k 

keN JQ<t 1 <...<t k <t 



J2 (m\o\m) (fi) 



where the first sum is taken over the number k of jumps, 
the integrals are taken over the jump times U, the remain- 
ing sums are taken over all possible types of jumps, and 
the \i/j(t)) is the unnormalized Monte Carlo wavefunction 
obtained from the initial wavefunction by the determin- 
istic non hermitian Hamiltonian evolution interrupted at 
times ti by the action of the jump operators of type a^. 
Here, for one body losses, the jump operators may be 
taken as C r = dy 1 / 2 r 1 / 2 -0(r), where r is any point on 
the grid of the lattice model (of unit cell volume dV). 
The jump associated to C r then describes the loss of a 
particle in point r. The non hermitian Hamiltonian is 
H cS = ff-|J r C$C r = H - ihTN/2, where N is the 
total number operator. 

The lossy ergodic model is based on a classical field 
model, where the state vector of the system is ap- 
proximated by a Fock state with N(t) particles in the 
mode cf>(r) linked to the classical field ip(r) by ip{r) = 
7V 1 / 2 (t)0(r). Then the action of C r on this Fock state 
simply pulls out a factor dV^-^T 1 / 2 ^^) in front of a Fock 
state with N(t) — 1 particles in the mode 4>. One may 
thus easily take the sum over the types of jumps in (|F1[) . 
oti corresponding to a loss event in rj, which produces 
factors equal to Y times the updated atom number after 
successive jumps. Also, in the lossy ergodic model, the 
condensate accumulated phase tp(t) is a classical quan- 
tity, evolving with the rate <p(t) = —fi mc [E(t),N(t)]/h. 

At zero temperature, one then simply has tp(t) = 
—gN(t)/hV so that, after a sequence of k jumps at times 



<p(t) 



fJ 



[Nh + (N- l)(t 2 -*i) + 



nv 

+(N-k)(t-t k 



9 

hV 



(N - fc)f + 53 



(F2) 
(F3) 



where N is here the initial atom number N(0). Similarly 
the squared norm of the Monte Carlo wavefunction after 
that sequence of jumps is 



(il>(t)\il){t)) = T k N(N-l)...{N-k + l) 



x e 



-r[JVti+...+(JV-fe)(t-t k )] 



Nl 



3 -r[(jv-*)t+£? =1 u 



(F4) 
(F5) 



(N — k)\ 

The expectation value of tp n (t), n positive integer, is thus 

N . k 

(<p n (t)) = V C k N Y k e- T ( N -W / dh... dt kl p n (t) TT , 

k~ -Am* ti 

(F6) 

where = N\/[kl(N — fc)!] is the binomial coefficient 
and where we used the fact that the integrand was a 
symmetric function of the times ti to extend the integra- 
tion domain to [0,t] h after division by k\. After lengthy 
calculations, and using the values of the binomial sums 



-rti 



^C k N e xk k n = d n x {e x + l) N , 



k=0 



(F7) 



we obtain the zero temperature results of the model: 
gN 



hvv 



(1 - e-™) 



Var^(<) = yL:) N(l-2Tte 



n 



(F8) 
2rt ) .(F9) 



This gives (|65)) . 

Next, we solve the lossy ergodic model to first order 
in r, at a non-zero temperature. To this order, one 
can restrict to the contributions of the zero-jump and 
of the single-jump trajectories. Calling /i the initial mi- 
crocanonical chemical potential, a function of the initial 
(random) energy E and (fixed) atom number N, and call- 
ing fi + Sfi the value of the chemical potential after the 
first jump, we have —tup(t) = fit for the zero-jump trajec- 
tory and —fup(t) — fiti + (fj,+8fi)(t — ti) — fit + 6 fi(t — 1\) 
for the single-jump trajectory with a jump at time t\. 
Thus 



<[-M*)] n > = M n )e 
xTNe~ r ( N 



-TNt 



dh {( 



-TNt! 



w-tj^t + s^t-tiW^ + oir 2 ), 

(F10) 
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where (. . .) in the right hand side stands for the expec- 
tation value over the initial system energy E. To first 
order in T, the exponential factors in the integral may be 
replaced by unity. Performing the integral over t\ gives 



{-fhp(t)) 
{[-tup{t)f) 



(p)t+(5p) l -TNt 2 + 0{T 2 ) (Fll) 



(M 2 )t 2 



- {{5p) 2 ) l -TNt 3 



(p5p)TNt 3 + o(r 2 ). 



(F12) 



This leads to (|63p . Note that the final result here is valid 
for Tt -C 1, even if our derivation seems to request the 
stronger condition TNt <C 1. 

Explicit expressions may be obtained from 
(inH]),!!!!]),®,©,®, and in the thermodynamic 
limit, where in particular one may approximate f{N — 1) 
by f (N) - df(N)/dN. Setting 



M ^ V 2m V ) 



(F13) 



m) = §j( N) + V s2{N) + lv s{N) > (F14) 



perature limit: 



S(N) 



M w VK 3 
1 



fs 2 (N) 



(F19) 
(F20) 

< S(JV)(F21) 



« 5(JV) (F22) 



1 



S(iV) 



N NknT 



S{N)H{N)Y&v{E - E ) ps — 

P 

{E-E )E(N) ps — . 



(F23) 
(F24) 
(F25) 



Since K 3 / p is of the order of the non condensed fraction 
(SN)/N, supposed to be <C 1 here, we recover 



APPENDIX G: QUANTUM SINGLE MODE 
MODEL WITH ONE BODY LOSSES 



and noting that S(N) = 0(V°) and S(iV) = 0(1/V) in 
the thermodynamic limit, we have 



H = f?[N + (E - E )S(N)} 



SE-SE, = - E ~ E ° ■ ^ • r,nr-i 



(F15) 

- (p - /xo) + 0(^- A )(F16) 



5/i 



V 



N 

[l + (E-EoMN)] 



+ o(v- 2 ) 



(F17) 



where 5E is the energy change after the first jump and 
//q = gN/V is the zero temperature (classical field) chem- 
ical potential. Taking the expectation value in (|63|) over 
the initial system energy E gives 

Yaxhp{t) ~ S a (JV)* 2 Var(^-Eb) 

+ IrTVt 3 2 [1 - 35(iV)E(iV)Var (£ - B ) 
+ 2(E - Eq)Z(N) + (E- E q ) 2 J: 2 (N)] . (F18) 

Here one simply has Var (E—Eq) = Var E since the initial 
particle number is fixed so that the ground state energy 
Eq does not fluctuate. For a classical field model in the 
canonical ensemble, (E — Eq) = M.k B T and Var (E — 
E )=M(k B T) 2 . 

In the limit k B T ^> Ng/V, which is natural for a clas- 
sical field model, the above expression for Var fup(t) may 
be greatly simplified. Taking a momentum cut-off K such 
that h 2 K 2 /2m = k B T, and ignoring numerical factors, 
we obtain in the thermodynamical limit and high tem- 



We show here that (|55|) . obtained at zero temperature 
within a classical field model, extends to the quantum 
case of a pure condensate with a large atom number and 
in an initial number state with N particles. 

The master equation for the single mode quantum 
model density operator p with one body losses is 



d P 1 rr> -l 

dt = m [H > p] 



Taop&l 



2 {°i & o, />} 



(Gl) 



where do annihilates a particle in the condensate mode, 
and H = ga£dl/{2V). A useful consequence is that the 
mean value of a not explicitly time dependent operator 
O evolves as 



dt 



(6) 



dt 



Tr[6p] 



,1 

r 



[o,H]) 

[al[6,a ] + [ald]a ). (G2) 



Neglecting the possibility that the condensate mode is 
empty, we use the modulus-phase representation clq = 
e^N^ 2 where the phase operator 9 and the number 
operator Nq — a^do obey the commutation relation 
[6, No] = — i- I n Heisenberg picture, the incremental 
evolution of the phase operator during an infinitesimal 
time step dt involves, in addition to the usual commuta- 
tor with the Hamiltonian H, a deterministic term A and 
a quantum stochastic term dB scaling as dt 1 / 2 due to the 
losses (H: 



dO 



dt 
ih 



[0,H]+ Adt + dB. 



(G3) 
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Applying (|G2|) to O = 8 gives 

1 = 0. 

Applying $g2§ to 6 = 9 2 gives 

(dB 2 ) = Tdt([al,6][6,a ] 



'Wo' 



(G4) 



(G5) 



In the large occupation number limit, we may thus ne- 
glect dB and take 



dt 



6~[d,H]/ih = -g(N -l/2)/hV. (G6) 



This justifies the assumption in the classical field model 
that the condensate phase is not affected by a jump. In 
the quantum model, the variance of the condensate ac- 
cumulated phase ip(t) — L dr[d§ (r) / dr] is thus 



dr / dr' 



(N (t)N (t')) 



{Mr)) (Mr')) ■ (G7) 



To calculate the one time averages of Nq and Aq we 
use (IGS with 6 = N and 6 



A 2 : 



d(N )/dt 
d(N 2 )/dt 



-T(N ) (G8) 
-2r(A 2 )+r(A ) (G9) 



that are straightforward to integrate with initial condi- 
tions (Mt = 0)) = N and (N 2 {t = 0)) = A 2 . 



To calculate the two time averages, we can restrict to 
r > t' by hermitian conjugation. Then we use the quan- 
tum regression theorem: setting <t(t') = A (0) j 5(r'), the 
"density operator" <t(t) evolves at later times r > r' with 
the same master equation as p, and 

(Mr)Mr')) = Tr[A> (0)a(r)] (G10) 

for t > t' . As a consequence 

d 



dr 



(N (t)N (t')) = -T{N (t)N (t')) (Gil) 



for t > t' , which is straightforward to integrate with 
the initial condition at r = r', (Nq(t')). We obtain for 
r > t' > 0, and for an initial number state with N 
particles: 

(Mr)) = Ne- TT (G12) 
<7V 2 (r)) = N 2 e- 2FT + Ner rT (l-e- r \pl3) 
(Mr)Mr')) = e- r ^(N 2 (r')). (G14) 



Mapping the double integral in (|G7|) to the integration 
domain < r' < r leads to 



Var 



1 - 2Re~ r * - e" 2rt l 



(G15) 



which coincides with the zero temperature classical field 
model result (1631). 
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